#include "mloop.h"
#include "grid.h"

void grid_mapping ()
{
	int i, j, k, l, n;
	float cx, cy, cz;
	float minx, maxx, miny, maxy, minz, maxz;
	int x1, x2, y1, y2, z1, z2;
	int grid1D[MAX_GRID_LINE][4];
	int *p;

	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
				grid.cell[i][j][k].number = 0;

	n = 0;
	for (i = 0; i < atomnum; i++)
	{
		minx = environatom[i].coor[0] - R_MAX - environatom[i].radius - grid.center[0];
		maxx = environatom[i].coor[0] + R_MAX + environatom[i].radius - grid.center[0];
		miny = environatom[i].coor[1] - R_MAX - environatom[i].radius - grid.center[1];
		maxy = environatom[i].coor[1] + R_MAX + environatom[i].radius - grid.center[1];
		minz = environatom[i].coor[2] - R_MAX - environatom[i].radius - grid.center[2];
		maxz = environatom[i].coor[2] + R_MAX + environatom[i].radius - grid.center[2];

		minx = GRIDMAX_X * 0.5 + minx / GRID_SIZE;
		maxx = GRIDMAX_X * 0.5 + maxx / GRID_SIZE;
		miny = GRIDMAX_Y * 0.5 + miny / GRID_SIZE;
		maxy = GRIDMAX_Y * 0.5 + maxy / GRID_SIZE;
		minz = GRIDMAX_Z * 0.5 + minz / GRID_SIZE;
		maxz = GRIDMAX_Z * 0.5 + maxz / GRID_SIZE;

		if (maxx < 0 || maxy < 0 || maxz < 0)
			continue;
		x1 = (int) minx - (minx > 0 ? 0 : 1);
		y1 = (int) miny - (miny > 0 ? 0 : 1);
		z1 = (int) minz - (minz > 0 ? 0 : 1);
		x2 = (int) maxx + 1;
		y2 = (int) maxy + 1;
		z2 = (int) maxz + 1;

		x1 = x1 < 0 ? 0 : x1;
		y1 = y1 < 0 ? 0 : y1;
		z1 = z1 < 0 ? 0 : z1;
		x2 = x2 > GRIDMAX_X - 1 ? GRIDMAX_X - 1 : x2;
		y2 = y2 > GRIDMAX_Y - 1 ? GRIDMAX_Y - 1 : y2;
		z2 = z2 > GRIDMAX_Z - 1 ? GRIDMAX_Z - 1 : z2;

		for (j = x1; j <= x2; j++)
		{
			minx = grid.center[0] + GRID_SIZE * (j - GRIDMAX_X * 0.5);
			maxx = minx + GRID_SIZE;
			cx = environatom[i].coor[0];
			if (cx > maxx)
				cx = cx - maxx;
			else if (cx < minx)
				cx = minx - cx;
			else
				cx = 0;
			for (k = y1; k <= y2; k++)
			{
				miny = grid.center[1] + GRID_SIZE * (k - GRIDMAX_Y * 0.5);
				maxy = miny + GRID_SIZE;
				cy = environatom[i].coor[1];
				if (cy > maxy)
					cy = cy - maxy;
				else if (cy < miny)
					cy = miny - cy;
				else
					cy = 0;
				for (l = z1; l <= z2; l++)
				{
					minz = grid.center[2] + GRID_SIZE * (l - GRIDMAX_Z * 0.5);
					maxz = minz + GRID_SIZE;
					cz = environatom[i].coor[2];
					if (cz  > maxz)
						cz = cz - maxz;
					else if (cz < minz)
						cz = minz - cz;
					else
						cz = 0;
					if (cx * cx + cy * cy + cz * cz < (R_MAX + environatom[i].radius) * (R_MAX + environatom[i].radius))
					{
						grid1D[n][0] = i;
						grid1D[n][1] = j;
						grid1D[n][2] = k;
						grid1D[n][3] = l;
						grid.cell[j][k][l].number++;
						n++;
						if (n == MAX_GRID_LINE)
						{
							printf ("MAX_GRID_LINE is too small\n");
							exit (-1);
						}
					}
				}
			}
		}
	}

	p = grid.grid1D;
	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
			{
				grid.cell[i][j][k].p = p;
				p += grid.cell[i][j][k].number;
			}

	for (i = 0; i < n; i++)
	{
		j = grid1D[i][1];
		k = grid1D[i][2];
		l = grid1D[i][3];
		*grid.cell[j][k][l].p = grid1D[i][0];
		grid.cell[j][k][l].p++;
	}

	p = grid.grid1D;
	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
			{
				grid.cell[i][j][k].p = p;
				p += grid.cell[i][j][k].number;
			}
}

void grid2_flag (float coor[3])
{
	int i, j, k, l;
	float cx, cy, cz;
	float minx, maxx, miny, maxy, minz, maxz;
	int x1, x2, y1, y2, z1, z2;

	minx = coor[0] - PACKING_CUTOFF - grid2.center[0];
	maxx = coor[0] + PACKING_CUTOFF - grid2.center[0];
	miny = coor[1] - PACKING_CUTOFF - grid2.center[1];
	maxy = coor[1] + PACKING_CUTOFF - grid2.center[1];
	minz = coor[2] - PACKING_CUTOFF - grid2.center[2];
	maxz = coor[2] + PACKING_CUTOFF - grid2.center[2];

	minx = GRIDMAX_X * 0.5 + minx / GRID_SIZE;
	maxx = GRIDMAX_X * 0.5 + maxx / GRID_SIZE;
	miny = GRIDMAX_Y * 0.5 + miny / GRID_SIZE;
	maxy = GRIDMAX_Y * 0.5 + maxy / GRID_SIZE;
	minz = GRIDMAX_Z * 0.5 + minz / GRID_SIZE;
	maxz = GRIDMAX_Z * 0.5 + maxz / GRID_SIZE;

	if (maxx < 0 || maxy < 0 || maxz < 0) return;
	x1 = (int) minx - (minx > 0 ? 0 : 1);
	y1 = (int) miny - (miny > 0 ? 0 : 1);
	z1 = (int) minz - (minz > 0 ? 0 : 1);
	x2 = (int) maxx + 1;
	y2 = (int) maxy + 1;
	z2 = (int) maxz + 1;

	x1 = x1 < 0 ? 0 : x1;
	y1 = y1 < 0 ? 0 : y1;
	z1 = z1 < 0 ? 0 : z1;
	x2 = x2 > GRIDMAX_X - 1 ? GRIDMAX_X - 1 : x2;
	y2 = y2 > GRIDMAX_Y - 1 ? GRIDMAX_Y - 1 : y2;
	z2 = z2 > GRIDMAX_Z - 1 ? GRIDMAX_Z - 1 : z2;

	for (j = x1; j <= x2; j++)
	{
		minx = grid2.center[0] + GRID_SIZE * (j - GRIDMAX_X * 0.5);
		maxx = minx + GRID_SIZE;
		cx = fmaxf (fabsf (coor[0] - minx), fabsf (coor[0] - maxx));
		for (k = y1; k <= y2; k++)
		{
			miny = grid2.center[1] + GRID_SIZE * (k - GRIDMAX_Y * 0.5);
			maxy = miny + GRID_SIZE;
			cy = fmaxf (fabsf (coor[1] - miny), fabsf (coor[1] - maxy));
			for (l = z1; l <= z2; l++)
			{
				if (grid2.cell2[j][k][l].flag == 1)
					continue;
				minz = grid2.center[2] + GRID_SIZE * (l - GRIDMAX_Z * 0.5);
				maxz = minz + GRID_SIZE;
				cz = fmaxf (fabsf (coor[2] - minz), fabsf (coor[2] - maxz));
				if (cx * cx + cy * cy + cz * cz < PACKING_CUTOFF * PACKING_CUTOFF)
				{
					grid2.cell2[j][k][l].flag = 1;
				}
			}
		}
	}
}

void grid2_mapping ()
{
	int i, j, k, l, n;
	float cx, cy, cz;
	float minx, maxx, miny, maxy, minz, maxz;
	int x1, x2, y1, y2, z1, z2;
	int grid1D[MAX_GRID2_LINE][4];
	int *p;

	grid2.center[0] = grid.center[0];
	grid2.center[1] = grid.center[1];
	grid2.center[2] = grid.center[2];

	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
				grid2.cell2[i][j][k].number = grid2.cell2[i][j][k].flag = 0;

	for (i = 0; i < atomnum; i++)
	{
		if (environatom[i].status & IS_MAINCHAIN_AND_CB)
			grid2_flag (environatom[i].coor);
	}

	n = 0;
	for (i = 0; i < atomnum; i++)
	{
		if (!(environatom[i].status & IS_MAINCHAIN_AND_CB))
			continue;
		minx = environatom[i].coor[0] - PACKING_CUTOFF - grid2.center[0];
		maxx = environatom[i].coor[0] + PACKING_CUTOFF - grid2.center[0];
		miny = environatom[i].coor[1] - PACKING_CUTOFF - grid2.center[1];
		maxy = environatom[i].coor[1] + PACKING_CUTOFF - grid2.center[1];
		minz = environatom[i].coor[2] - PACKING_CUTOFF - grid2.center[2];
		maxz = environatom[i].coor[2] + PACKING_CUTOFF - grid2.center[2];

		minx = GRIDMAX_X * 0.5 + minx / GRID_SIZE;
		maxx = GRIDMAX_X * 0.5 + maxx / GRID_SIZE;
		miny = GRIDMAX_Y * 0.5 + miny / GRID_SIZE;
		maxy = GRIDMAX_Y * 0.5 + maxy / GRID_SIZE;
		minz = GRIDMAX_Z * 0.5 + minz / GRID_SIZE;
		maxz = GRIDMAX_Z * 0.5 + maxz / GRID_SIZE;

		if (maxx < 0 || maxy < 0 || maxz < 0)
			continue;
		x1 = (int) minx - (minx > 0 ? 0 : 1);
		y1 = (int) miny - (miny > 0 ? 0 : 1);
		z1 = (int) minz - (minz > 0 ? 0 : 1);
		x2 = (int) maxx + 1;
		y2 = (int) maxy + 1;
		z2 = (int) maxz + 1;

		x1 = x1 < 0 ? 0 : x1;
		y1 = y1 < 0 ? 0 : y1;
		z1 = z1 < 0 ? 0 : z1;
		x2 = x2 > GRIDMAX_X - 1 ? GRIDMAX_X - 1 : x2;
		y2 = y2 > GRIDMAX_Y - 1 ? GRIDMAX_Y - 1 : y2;
		z2 = z2 > GRIDMAX_Z - 1 ? GRIDMAX_Z - 1 : z2;

		for (j = x1; j <= x2; j++)
		{
			minx = grid2.center[0] + GRID_SIZE * (j - GRIDMAX_X * 0.5);
			maxx = minx + GRID_SIZE;
			cx = environatom[i].coor[0];
			if (cx > maxx)
				cx = cx - maxx;
			else if (cx < minx)
				cx = minx - cx;
			else
				cx = 0;
			for (k = y1; k <= y2; k++)
			{
				miny = grid2.center[1] + GRID_SIZE * (k - GRIDMAX_Y * 0.5);
				maxy = miny + GRID_SIZE;
				cy = environatom[i].coor[1];
				if (cy > maxy)
					cy = cy - maxy;
				else if (cy < miny)
					cy = miny - cy;
				else
					cy = 0;
				for (l = z1; l <= z2; l++)
				{
					if (grid2.cell2[j][k][l].flag == 1)
						continue;
					minz = grid2.center[2] + GRID_SIZE * (l - GRIDMAX_Z * 0.5);
					maxz = minz + GRID_SIZE;
					cz = environatom[i].coor[2];
					if (cz  > maxz)
						cz = cz - maxz;
					else if (cz < minz)
						cz = minz - cz;
					else
						cz = 0;
					if (cx * cx + cy * cy + cz * cz < PACKING_CUTOFF * PACKING_CUTOFF)
					{
						grid1D[n][0] = i;
						grid1D[n][1] = j;
						grid1D[n][2] = k;
						grid1D[n][3] = l;
						grid2.cell2[j][k][l].number++;
						n++;
						if (n == MAX_GRID2_LINE)
						{
							printf ("MAX_GRID_LINE is too small\n");
							exit (-1);
						}
					}
				}
			}
		}
	}

	p = grid2.grid1D;
	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
			{
				grid2.cell2[i][j][k].p = p;
				p += grid2.cell2[i][j][k].number;
			}

	for (i = 0; i < n; i++)
	{
		j = grid1D[i][1];
		k = grid1D[i][2];
		l = grid1D[i][3];
		*grid2.cell2[j][k][l].p = grid1D[i][0];
		grid2.cell2[j][k][l].p++;
	}

	p = grid2.grid1D;
	for (i = 0; i < GRIDMAX_X; i++)
		for (j = 0; j < GRIDMAX_Y; j++)
			for (k = 0; k < GRIDMAX_Z; k++)
			{
				grid2.cell2[i][j][k].p = p;
				p += grid2.cell2[i][j][k].number;
			}
}
